Principles of Numerical weather prediction and modeling 


There are three types of weather forecasts. They are short range, medium range 
and long range weather forecasting. 


Short range = few hours to 48 hours or sometimes 72 
hours 

Medium range = 3 days to about 3 weeks 

Long range = a few seasons to a year 


There are three weather forecasting methods. They are i) synoptic weather 
forecasting ii) Numerical methods and iii) statistical methods. 

Synoptic forecasting: 

Synoptic means the observation of different weather elements refers to a specific 
time of observation. From the careful study of weather charts over many years, 
certain empirical rules can be formulated. These rules will help to estimate the rate 
and direction of movement of weather systems like a trough line, hurricane or 
thunderstorm. 

Numerical methods: 

Numerical method is based on the fact that gases of the atmosphere follow a 
number of physical principles. If the current conditions are known, these physical 
laws help to forecast the future condition. 

A series of mathematical equations are used to develop theoretical models of 
the general circulation of the atmosphere. These equations are used to specify 
changes in the atmosphere as time changes using the weather elements like wind, 
temperature, humidity, clouds, rain, snow, evaporation etc. 

The atmosphere is divided into 6 or 11 layers for this purpose. As the initial 
conditions ( initial state of the atmosphere) are known through observations using 
surface and upper air data, it is possible to predict using numerical models for 
future conditions. 

Valuable contribution in the development of numerical methods was made by 
Prof.J.Charney in 1948 and Obukov. 

As there is large number of weather variables and only a less number of 
equations are available, some weather variables are ignored on the assumption that 
they do not change with time. As large numbers of calculations are involved even 
with simplified models, a high speed super computer is necessary. 

For carrying out numerical models the following considerations are followed: 
The equations are simplified 
The equations are designed in such a way that they guarantee the conservation of 
air mass, momentum, water vapor and total energy for all time. 

The region to be forecasted is divided into number of nodal points in a rectangular 

shape. 

The equations are solved at each nodal point for a short period of 10 minutes. 

By repetitive calculations for every next 10 minutes forecast is obtained for 24, 48 

or 72 hours. 

L.F.Richardson in 1922 attempted to solve the differential equations of momentum, 
energy, state and conservation using finite difference method with the help of desk 
calculator. 

After World War II and the development of digital computer, it was realized 
that Richardson’s scheme was not simple as the governing equations contains not 


only slow moving meteorologically important motions but also high speed sound 
and gravity waves. Though these waves are very weak, they amplify surprisingly in 
numerical methods and thus generate ‘noise’. 

So in 1948 J.G.Charney simplified the equations by introducing the scale 
analysis, geostrophic and hydrostatic assumptions. By these techniques sound and 
gravity waves are filtered. Then these equations are called quasi-geostrophic 
models. A special case if this model , called equivalent barotropic model, was used 
in 1950 which gave the first successful numerical forecast. 

A. Objective Analysis 

The automatic process of transforming of irregularly distributed observed data from 
different geographic locations at various times to numerical values at regularly 
spaced grid points at fixed times is called objective analysis. In the earlier days of 
numerical weather prediction, the two-dimensional fitting of polynomial was applied 
to observational data in an area surrounding a grid point at which the value of 
analysis is required. Polynomial fitting is suitable for processing relatively dense and 
redundant data, but it often gives unreasonable results in data of sparse areas. 
More recent procedures of objective analysis practiced at most of the operational 
forecasting centers consist of the steps schematically described in Figure below. 
Most observed data are collected worldwide every 6 hours, say 00, 06,12 etc. UTC 
(Coordinated Universal Time), and are blended into forecast values appropriate to 
the map time by the process referred to as data assimilation. The basic process of 
data assimilation, indicated by the box enclosed by dashed lines in Figure below, 
consists of two steps: One is objective analysis and the other is called initialization. 
The step of objective analysis is carried out by statistical interpolation by taking into 
account all of the available observations plus other prior information. Short-term 
forecasts from the previous analysis cycle are used as prior information and are 
referred to as background fields. 

B. Concept of Initialization 

The solutions of primitive equation models correspond to two distinct types of 
motion. One type has low frequency. 

Its motion is quasi-geostrophic and meteorologically dominant. The other 
corresponds to high-frequency gravity-inertia modes. The amplitude of the latter 
type of motion is small in the atmosphere. Hence, it is important to ensure that the 
amplitudes of high-frequency motions are small initially and remain small during the 
time integration when the primitive equations are solved as an initial value problem. 
The process of adjusting the input data to ensure this dynamical balance is called 
initialization. 

Solution to the initialization problem was central to the successful transition in 
forecast practice from the use of 

quasi-geostrophic models to primitive equation systems during the 1970s. Since 
gravity-inertial motions are filtered 

out in quasi-geostrophic models, no special procedure was necessary and the 
objectively analyzed data could be used immediately as the input data to quasi- 
geostrophic models. Actually, there is an interesting history that led to solution of 
the initialization problem that is equivalent to the quest of understanding what the 
primitive equation forecast system really is. A promising break came with the 
advent of so-called nonlinear normal mode initialization (NNMI) during the latter 
part of the 1970s. 
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Numerical solution of equations 


The five prognostic equations used for forecasting are: 
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The first three equations are equation of motion and the other two are equations of 
continuity and law of conservation of energy for an adiabatic process. There are five 
dependent variables u,v.w,p and op and five equations. So they can be solved. 


Finite difference scheme to solve the equations 


Using the Taylor’s equation we can write 
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Taking first two terms in the right hand side, we can write 
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The finite difference grid is as below: 
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For example, to get ay at middle point ‘O’ we can write: 
Op _ Pr-p, 
dy 2d 


Similar expressions can be written for all other five variables. Thus incorporating all 
these values at each grid point, future time value is known. 


Filtering of Sound and Gravity waves 


Sound waves are longitudinal waves which are generated because of 
alternating adiabatic compression and expansion of the medium. This arises 
because of pressure gradient. In the case of hydrostatic atmosphere, the vertical 
pressure gradient can not be influenced by adiabatic compression. Hence the 
vertically propagated sound waves are filtered under the hydrostatic approximation. 


However, horizontally propagating sound waves are eliminated by the assumption 
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After the sound waves are filtered, the prediction equations, such as momentum 
equation, continuity equation and hydrostatic thermodynamic equation are written 
as: 
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equations along with the hydrostatic and equation of state are called primitive 
equations. They contain the gravity waves as solutions. 

The gravity waves are waves in which buoyancy acts as the restoring force 
on air parcels displaced vertically. A disturbance is propagated because of alternate 
horizontal convergence and divergence in individual columns of fluid. A divergent 
horizontal velocity field which changes in time is essential for the propagation of 
gravity waves. Neglecting the local rate of change of the horizontal divergence in 
computing the balance between mass and velocity fields is sufficient to filter out the 
time dependent gravity waves. As the prediction equations of 1 to 3 do not contain 
the local derivatives V. V explicitly, we use the vorticity and divergence equations 
instead of horizontal momentum equations. 

The vorticity and divergence equations are as below: 
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The divergence equation is 
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By operating V.() on this equation, the divergence equation can be written as : 
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If the Left hand side of equation (5) is set to zero, all solutions corresponding to 
time dependent gravity waves can be eliminated. But the equation (5) becomes an 
intractable diagnostic relationship between , VA\q . Hence we do the scale analysis 


for systematic elimination of unimportant terms. Usually neglecting the second 
order terms is sufficient to filter the gravity waves. 

In order to exhibit the connection between scaling and filtering, it is 
convenient to partition the wind as 


V=V,+V, 
Where V, is the non divergent part and V. is the divergent part of the wind. 
If the velocity field is two dimensional, the non-divergent part can be 
expressed in terms of stream function as 
Vj=KxVyp 


As V.V,=0, Vx V,=0 and 


The horizontal divergence is small compared to the vorticity (by scale analysis) 
for synoptic scale systems in mid latitudes. That is V is quasi- nondivergent 


Vy>lVe 


Thus to a first approximation we can replace V by V, everywhere in equations 4 
and 5except in terms involving horizontal divergence. 
So the resulting vorticity and divergence equations after scaling are 
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In the above equations gravity wave terms were filtered. From equation (10), we 
can interpret as 
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This means the geopotential field on a constant pressure chart is approximately 
proportional to the stream function field. Equation 10 also says that to a first 
approximation, the vorticity is the vorticity of geostrophic wind calculated using the 
constant mid latitude value of the coriolis parameter. 


The geostrophic vorticity equation and the hydrostatic thermodynamic 
equation can be written as: 
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By eliminating w from these equations, the quasi-geostrophic potential vorticity 
equation as: 
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We have assumed that a is a function of pressure only. This relation states that the 
geostrophic potential vorticity is conserved following the non-divergent wind in 
pressure coordinates. 
In practice solving equation (13) is extremely difficult and so one has to go for 
numerical schemes. Models have been developed which simplify the numerical 
structure by referring to one, two or few layers. 
One Parameter models 
There are some simplest models with only one cluster level in the vertical. 
The barotropic model: 


In this model one assumes that there exists a level of non divergence in the 
atmosphere. Usually this level is assumed as 500mb level. 


Thus JV. vee =0 at 500 mb level 


Hence equation (9) becomes =| Va y\=- V,- Z| a y+ £| betas iaeuineadas (14) 


This equation allows to compute the evolution of the fluid at the level of non 
divergence only. 


b) The equivalent barotropic model: 


Same type of variation of wind with height is taken in equivalent barotropic 
model by averaging in the vertical. But one level is only considered in the 
model . In this model the nondivergent wind is expressed as: 
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The angle bracket denotes the vertical average and A(p) is a weight function. 
The vorticity equation for this model is written at the stairred level where 
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CA w are also expressed as, 
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To the extent that equation (15) is correct, the results from the solution of 
vorticity equation can be utilized to deduce the evolution of yw at other levels. 


Two parameter baroclinic model 


The one parameter models do not allow temperature advection. But the thermal 
advection processes are essential to the development of synoptic systems. The 
barotropic models are quite effective in predicting the evolution of mid tropospheric 
flow for periods upto a couple of days. For consistently reliable forecast one should 
predict the development of new systems also. Hence there is a need of baroclinic 
models. The simplest baroclinic model is a two layer model. The atmosphere is 
divided into four discrete layers and the disposition of variable is given in the 
following schematic diagram. 


py = 0 hPa 


pi = 250 hPa 


po = 500 hPa 
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ps4 = 1000 hPa MU 


Here ‘I’ is level and p is pressure surface. 


For a flat ground the lower boundary condition is w, = 0. Applying the vorticity 


equation (11) at levels 1 and 3 and utilizing the finite difference technique, we 
can write: 
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The vorticity equations at the two levels are: 
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Now equations 16, 17 and 18constitute a close set of prediction equations. 


First eliminate w, between 16 and 17 by adding both equations to get 
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This equation governs the barotropic part of the flow. 
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Next subtract ( 17) from (16 )and add the results to ae times (18) to get 
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Here yu is anon dimensional parameter. Equation (19) states that the 
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local rate of change of vertically averaged vorticity is equal to the average of 
vorticity advections at the two levels. Equation (20) gives the local rate of 
change of thickness of 250-750 mb layer. 


The time derivatives between equations 16, 17 and 18 can be eliminated to get 
the omega equation as 
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Since equations 19, 20 and 21 are derived from the quasi-geostrophic system of 
equations, in the two parameter model, the vorticity changes are geostrophic and 
the temperature changes are hydrostatic. Due to these constraints, the omega field 
is to be adjusted at every instant. 


Primitive Equation models 


The filtered models discussed above, involves various approximations of the quasi- 
geostrophic system. So there are some limitations on the accuracy of forecasting. In 
low latitudes, the stream function yw can not be obtained from the geostrophic 
vorticity field and can be obtained from more complicated balance equation which 
is a nonlinear equation. It is difficult and time consuming to solve the balance 
equation at each time step. 

The following balance equation is derived from the divergence equation by 


putting initially < V.V\=0 and 7. V=0. 
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Therefore, it is worthwhile to investigate the possibility of using less restrictive 
modeling assumptions. 
The scale analysis indicated that the hydrostatic approximation is far more accurate 
than other filtering approximations. Hence we continue to assume that the motions 
are hydrostatic. The resulting system of equations is isobaric coordinates consists of 
three prognostic equations (x and y components of momentum equations and 
thermodynamic equation) and three diagnostic equations (the continuity, 
hydrostatic ad the equation of state). 

These six equations constitute a closed set with the dependent variables 
u,v,w,@, 8, @ which are referred to as primitive equations. 

Let us write the basic primitive equations in isobaric coordinate system as: 
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The transformation from p to o coordinates can be derived in a similar manner 
of transformation from Z to P coordinates. 


Isobaric and sigma coordinates 


In isobaric coordinates we put WP |=—P)gw|Z,|. We assume that the height 


of the ground Z, is coincident with the pressure surface po ( = 1000mb). However, 
this is not strictly true, as pressure surface will not coincide with the ground even 
when the ground is level. So there is a need of o system of coordinate is used in the 


vertical in numerical modeling such that aa where p; denotes the pressure at the 


Ss 
surface. In sigma coordinate system the lower boundary condition is applied at 
o=1. 
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The vertical velocity in o coordinates is 0 = ae will always be zero at the ground 


as o=1 is constant even in the presence of sloping terrain. Hence in o coordinates 
the lower boundary condition is 0 =Oat o=1. 


0 
Using the definition of , one can write: V,()= Ale Vaal 


So the system equations in o coordinates can be written as: 
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The set of six equations in (23) contains seven dependent variables u, v, 0, 6, Tq 
and p,. So another equation is necessary to solve these variables. This equation is 
surface pressure tendency equation. It can be obtained by integrating the continuity 
equation vertically by using the boundary conditions 0 =O at o=OA1. 
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This equation states that the rate of increase of the surface pressure at a given 
point is equal to the mass convergence in a unit cross section of column above that 
point. 


The surface pressure tendency equation is 


Two layer primitive equation model (PE Model 


In the primitive equation model the static stability parameter o is a function 
of space and time whereas in earlier models, o(p) is usually a specified parameter. 


It can be shown that if s is specified as a function of pressure alone, that is 

o=o(p), it is not possible to formulate a PE model which satisfies the requirements 
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of energy conservation. Hence it is necessary to compute Op explicitly at each 


time step. Hence temperature must be predicted at a minimum of two levels rather 
than at one level as in the case of Z-level quasi-geostrophic model. 
The schematic diagram of Z-level PE model with variables is as below: 
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The momentum equation and the thermodynamic equations are applied at levels 
1 and 3. The vertical differencing of these equations requires knowledge of 
sigma vertical motion at level 2. GO, is obtained diagnostically through the use 
of continuity equation as follows: 


JO 0O,—Oo : 
Ps 
($2) al a 
2 
ag\ _G-G%_ a. 
Sal —— a 


From the continuity equation in (28) we can write: 
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Adding (24) and (25) we can write: 
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Subtracting (25) from (24) we can get: 
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In asimilar manner @, and @, can be obtained diagnostically by vertical 
integration of hydrostatic relation in equation (23). 


In short the following steps may be followed for the forecast using PE model. 


1. Suitable finite difference analogues of the momentum and thermodynamic 
equations at levels 1 and 3 and the surface pressure tendency equation are 
written. 


2. These prognostic equations are utilized to obtain tendencies of V,, V;,0,,03 
and P.. 

3. Extrapolating the tendencies ahead using a suitable time differene scheme. 

4. The new variables of V,, V3,0,,0;, and p, are used to determine 6;, @, and 
@; diagnostically. 

5. Repeat steps 2 to 4 until the desired forecast time is reached. 


It must be noted that the PE model contains the mechanism for both gravity wave 
and horizontal acoustic wave propagation. It is necessary that the time increment 


co bes. where C is horizontally 


be kept small enough to satisfy the condition: 


propagating sound wave which is the fastest moving wave solution, d is horizontal 
resolution and At is time increment. Therefore the time steps in PE models are 
considerably small ( =10 minutes). 


Initialization Scheme 


Initialization scheme involves the use of accurate initial data. One reason for the 
failure of Richardson is due to lack of accurate initial data. Though the data 
coverage is vastly improved since then actual measured wind and pressure data is 
still lacking as most of the data is derived from remote sensing satellite data. 


To introduce the initialization method, consider the relative magnitudes of 
the various terms in the momentum equation in isobaric coordinates: 
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For synoptic scale motions in mid latitudes, we have seen that the wind and 
pressure fields are in approximate geostrophic balance. Thus the acceleration 
following the motion is measured by the small difference between the two nearly 
equal terms Coriolis and geopotential (—fK x VA V@i. Although @can be 
determined with good accuracy using observed data, observed winds are often 10 
to 20% in error. Hence if wind is used as initial data, the Coriolis force works out to 
10 to 20% in error. Since the acceleration is normally 10% of the Coriolis force, the 
acceleration comes to 100% in error (10x10). Such spurious acceleration leads to 
poor estimates of initial pressure and velocity tendencies. This further produces 
large amplitude gravity wave oscillations. 


To avoid this problem, let us derive the balanced equation as follows.In the 


divergence equation at t = 0, put c+ V.V=Oand V.V=0 
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Since we have assumed that the initial velocity is non divergent, we can replace V 
by w by letting 


V,, = Kx Vw then we get 
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This equation is called balance equation. Equation (26) is linear in @ and non linear 
in yw. In midlatitudes @is known and henc e wy is solved with the help of equation 26. 


The whas two solutions. One solution corresponds to the gradient wind balance and 


the other one is anomalous. One y is known, the initial velocity field V,, can be 
computed. This velocity is in gradient wind balance with the geopotential field. The 
initial local accelerations are thus very small. Since the initial rate of change of 
divergence is zero, the large amplitude gravity wave oscillations are not excited by 


the initial conditions. 


